Mean field dynamics of superfiuid-insulator phase transition in a gas of ultra cold 

atoms 



in 

o 
o 

C 
(N 



5-H 



Ctf 



O 

m 
> 

00 

\o 
o 
^J- 

o 

-I— > 



i 

c 
o 
o 



X 
S3 



Jakub Zakrzewski 

Instytut Fizyki im. Mariana Smoluchowskiego and Mark Kac Complex Systems Research Center, 

Umwersytet Jagielloriski, ulica Reymonta 4, PL-30-059 Krakow, Poland 

(Dated: 2nd February 2008) 

A large scale dynamical simulation of the superfluid to Mott insulator transition in the gas of 
ultra cold atoms placed in an optical lattice is performed using the time dependent Gutzwiller mean 
field approach. This approximate treatment allows us to take into account most of the details of the 
recent experiment [Nature 415, 39 (2002)] where by changing the depth of the lattice potential an 
adiabatic transition from a superfluid to a Mott insulator state has been reported. Our simulations 
reveal a significant excitation of the system with a transition to insulator in restricted regions of the 
trap only. 
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I. INTRODUCTION 

A theoretical suggestion [J of a possibility to realize 
one of the standard models for interacting particles - the 
Bose-Hubbard (BH) model |3, Q] in a cold gas placed in 
an optical lattice has been followed soon by a seminal 
experiment |j|. The reported realization of a quantum 
phase transition between superfluid (SF) and Mott insu- 
lator (MI) phases showed convincingly that it was pos- 
sible to control experimentally parameters of the model 
practically at will. This triggered several studies involv- 
ing Bose condensate [jj la, 0, IM l9l llCM as well as, more 
recently Fermi-Bose mixtures [fl], [13, H3 placed on the 
optical lattices (the reference list must be not complete 
bearing in mind that more than 70 papers with "optical 
lattice" in the title are listed in the cond-mat archive last 
year only). 

At the same time a number of groups [14 . Ha . IWL 1171 Il8| 
tried to understand the details of the very first experi- 
ment yj to check the underlying physics. To imagine the 
difficulty in modeling the experiment let us recall that 
it involves about 10 5 interacting atoms (bosons) placed 
in the harmonic trap and the three dimensional (3D) 
lattice potential. Such a system is well described by a 
Bose-Hubbard model with position dependent chemical 
potential |l| . Even finding the ground state of the sys- 
tem for that number of particles and 65 x 65 x 65 lat- 
tice sites is a formidable task. State of the art quantum 
Monte Carlo (QMC) E 111 calculations aimed at 
the ground state properties include up to 16 sites in 3D 
|14J . more sites may be included in one (ID) or two (2D) 
dimensional models [ljall3- These studies, while inter- 
esting on their own, can shed little light on the dynamics 
of the system when its parameters are varied. Except for 
special exactly solvable models, the efficient simulation 
of time-dependent properties of interacting many-body 
system remains an open problem although recently quite 
a progress has been obtained for ID systems [2(J, |21| . 

It seems, therefore, that the only reasonable and 
tractable way of analyzing the dynamics of the discussed 
experiment is using approximate methods. To this end 



we shall use an approach based on the time dependent 
variational principle with Gutzwiller ansatz. That will 
allows us to model the details of the 3D experiment Q| . 
The prize for it is similar to that paid in other approx- 
imate treatments - one may always question the extend 
to which the approximations allow to describe the prop- 
erties of the systems studied. We hope to convince the 
reader that the numerical results are at least mutually 
consistent and thus may provide considerable insight into 
the dynamics of the experiment. 

The discussion of the dynamics is postponed to Sec- 
tion III since we discuss in the next Section the static 
mean field solutions for the ground state for experimen- 
tal parameters. Here a comparison with available exact 
QMC results is possible at least. This shall give us some 
confidence about the applicability of the mean field ap- 
proach yielding, at the same time, the initial state for the 
dynamics studied later. 



II. STATIC MEAN FIELD FOR THE 
BOSE-HUBBARD MODEL 

The Bose-Hubbard Hamiltonian describing the system 
takes the form ll 



h = -j J2 <4c 

<i,3> 



U 
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where Hi = a\di is an occupation number operator at site 
i (with a, being the corresponding annihilation bosonic 
operator), U the interaction energy, J the tunneling coef- 
ficient and Wi the energy offset at site i. J2<% 7 > denotes 
a sum over nearest neighbors. Both J and U are func- 
tions of the lattice potential and may be easily expressed 
in terms of integrals of the Wannier functions of the low- 
est energy band for cold atoms implementation of the 
model Q. 

Consider first the standard homogeneous situation in 
which all Wi's are equal. The last term in JIJ becomes 
proportional to the (conserved) number of bosons and 



may be dropped. The only remaining parameter of the 
model is the ratio of U/J. When tunneling dominates 
the system in its ground state is superfluid while in the 
opposite case it becomes the Mott insulator. The bor- 
derline between the two phases depends on the chemical 
potential /i. The Mott insulator state is incompressible 
and is characterized by an integer mean occupation of 
sites. In effect starting from the superfluid at small U/J 
and a non integer ratio of N/M (N denotes the number 
of atoms while M the number of sites) and increasing 
U/J the ground state remains superfluid up to highest 
values of U/J at fixed boson density. On the other hand 
the range of \x values corresponding to the commensu- 
rate filling increases with U/J. In effect, the separation 
line between a MI and a SF forms characteristic lobes 

For a detailed discussion of the BH model see [2], |16| . 
As we are interested in the mean field approximation, 
let us just quote Zwerger [l(j saying "In two and three- 
dimensional lattices, the critical value for the transition 
from a MI to a SF is reasonably well described by a mean- 
field approximation" . In one dimension the mean field 
approximation is much worse [16j |. 

In the presence of the additional potential, e.g. the 
harmonic trap, local energies Wi depend on the sites lo- 
cation. Then the effective chemical potential at each site 
becomes \x% = \x—Wi- As pointed out already in [lj this 
will lead for large U/J to a shell like structure with MI 
phases with different integer occupations (highest in the 
middle assuming attractive binding additional potential) 
separated by SF regions. This picture has been nicely 
confirmed in quantum Monte Carlo calculations both in 
ID and in 3D Q. 

The latter exact results are of particular interest for 
us since they allow for a comparison with the mean field 
approximation. In [14J a 3D 16 x 16 x 16 lattice is con- 
sidered with different values of U and J parameters as 
well as the harmonic trap. To find the mean field ground 
state we minimize 



< E>=<G\H -fiN\G>, 
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where N = J^. 
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and \G > is the Gutzwillcr trial func- 
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The number of parameters /„ depends on the number 
of sites (here 16 3 ) as well as the maximal occupation at 
a given site n m . The average maximal occupation at the 
center of the trap is 2 fot the data considered in jlij . 
Therefore, it is sufficient to take n m = 7. That yields 
a minimization procedure over 32768 parameters. Such 
a number of parameters must lead to a spurious local 
minima, unless a very good estimate exists for the initial 
set of fn s (i.e. the initial \G >). Fortunately such a 
guess is quite obvious and is often termed a local mean 
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Figure 1: Mean field particle density distribution (on-site fill- 
ing factors) as a function of the distance from the center of the 
trap measured in the units of the lattice constant a. Filled 
circles correspond to numerical results, lines are drawn to 
guide the eye. The parameters of the BH model match those 
quoted in Fig. 1 of 114 . see the text for numerical values. 
Comparison with the latter results obtained within "exact" 
quantum Monte Carlo reveals the accuracy of the mean field 
approximation. 



field approximation. Namely at each site i one takes a 
solution for /„ 's corresponding to the homogeneous BH 
model with the effective chemical potential yn = /i — Wi. 
Provided Wi changes smoothly from site to site, such an 
approach should be an excellent approximation to a full 
mean field solution. And indeed it is, we have found 
for the data discussed below that the initial and final 
< E > [see Eq. (J2J] differ by at most 2%; the number 
of iterations of standard Numerical Recipes minimiza- 
tion packages [22( is slightly bigger than the number of 
parameters. 

The results obtained are presented in Fig. ^ in the 
same form as the corresponding plot in |l4| to make 
the comparison easier. The values of parameters cor- 
respond to those taken in fl4l . The notation used by us 
differs slightly from that of [JTj , inparticular the tunnel- 
ing constant J is denoted as t in [Tj] . For the purpose of 
this figure all parameters are expressed in units of J, i.e. 
J = 1. To make the comparison with notation used in 
|14J easier, the coefficient in the term proportional to the 
number of atoms at a given site in J2J) i.e. the difference 
of the on-site energy Wi and the chemical potential \i is 
expressed as Wi — \i = — Uq + U/2 



Kx t , where Xj is a 



position vector of site i with x| being the square of the 
distance of the z-th site from the center of the harmonic 
trap measured in units of the lattice constant. With such 
a definition the parameters used in Fig. ^ take the val- 
ues: panel (a): U = 24., U = -11.08, k = 0.19531; 
(b): U = 32., U = -28.08, k = 0.19531; (c): U = 80., 
U = -65.0, k = 0.97656; (a): U = 80., C/ = -90.0, 
K = 1.03062; (a): U = 80., U = -120.08, re = 2.00375; 
and (a): J7 = 80., U = -150.0, k = 1.75781. A com- 
parison of Fig. ^ with Fig. 1 of [lj] indicates that, as far 
as the average occupation at different sites is concerned, 
the mean field solution is in excellent agreement with the 
quantum Monte Carlo results. 

Instead of the occupation at different sites one may 
take a look at the momentum distribution, i.e. the quan- 
tity closely related to that measured in the experiment 
(see U and the discussion below). The momentum dis- 
tribution is given by |14| 



n k = |<£(k)| 2 5>' 
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(4) 
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where k is the wavevector, 0(k) a Fourier transform of 
the Wannier site-function. The latter yields a broad 
bell-shaped background and provides merely information 
about the lattice. The relevant information about the 
atoms is contained in the Fourier transform of < ajaj >. 
In the mean field approximation this factorizes for differ- 
ent i,j < ajaj >~< a\ >< a,j >. Such a factorization 
seems quite drastic and one may expect significant dif- 
ferences between the momentum distributions obtained 
from QMC and within the mean field approximation. 
It is really not so, however, for bosons in a harmonic 
trap as visualized in Fig. [21 for the mean field. That fig- 
ure should be compared with QMC results presented in 
Fig. 2 of [l4|. Observe that differences appear only for 
panel (d), the exact results yield significantly broader 
momentum distribution. As discussed in [14( almost no 
SF fraction is present in the QMC result corresponding to 
panel (d) . Then the factorization must affect strongly the 
momentum distribution since in a vast majority of sites 
< a, >= 0. Clearly, however, as long as some SF frac- 
tion is present in the system the mean field momentum 
distribution closely resembles the exact quantum results. 

This apparent quite close agreement between QMC re- 
sults and the mean field approximation forl6xl6xl6 
lattice and about 10 3 atoms is very encouraging in view 
of realistic experimental conditions 4]. Here both the 
external potential changes less rapidly (the size of the 
lattice is now 65 x 65 x 65) and number of atoms exceeds 
10 5 thus one may expect that the mean field approxima- 
tion works even better. 

While the test described above have been taken for 
some chosen (by authors of 14] ) values of U, J, as well 
as the trap frequency, to simulate the experiment we have 
to determine first the relevant range of parameters. From 
now on we shall measure the quantities of dimension of 
energy in the units of the recoil energy of 87 Rb atoms 
for light with a wavelength A = 2ir/K — 852nm, i.e. 



E r = h 2 K 2 /2M, where M is the mass of 87 Rb atoms. 
The depth of the optical lattice V changes from to 22E r . 
Finding Wannier functions for different values of V [2Jj 
we evaluate the corresponding U(V) and J(V) values. 
The energy offset at each site Wi has two components 
in the experiment |4|. One is the harmonic magnetic 
trap potential (time-independent), another is due to the 
Gaussian intensity profiles of lattice creating laser beams. 
The latter may be also approximated by a harmonic term 
4] the corresponding frequency is then dependent on V. 

To find the mean field ground state for different V val- 
ues and the number of atoms of the order of 10 5 one needs 
to solve a minimization problem over 2 x 10 6 parameters 
(with n m = 7 as before) which is hardly manageable. 
One may, however, use the symmetry of the problem (cu- 
bic lattice combined with spherically symmetric trap) to 
significantly reduce that number. Let i,j,k count the 
sites in X,y,z directions, respectively with each index 
taking the values from —32 to n s = 32 (yielding 65 sites 
in each direction). Due to the ground state symmetry it is 
enough to consider only the sites with 0<i<j<k<n s 
which reduces the number of minimized parameters to 
about 48 thousands. Needless to say we have checked 
on the smaller 16x16x16 problem that the symmetry 
reduced problem yields the same ground state as the full 
minimization. 

The results obtained are presented in Fig. [3] and are 
practically indistinguishable from the initial guess i.e. a 
wavefunction coming from local mean field approxima- 
tion discussed above. The chemical potential jjl has been 
adjusted (for each case) to have the average number of 
atoms N =< N >= £\ < rn > around 10 5 . This leads 
to more than two atoms (on average) per site in the cen- 
ter of the trap. To characterize whether the state is closer 
to being superfluid or Mott insulator we define the super- 
fluid factor 7s p — 1/N^2 < o» >< o| >• This factor is 
zero for pure MI state (when < <Zj >= as each node is 
in a Fock state) and reaches unity for Poissonian statis- 
tics at each node. While obviously it is not a "proper" 
order parameter for the phase transition, it seems to be 
convenient for characterizing the states obtained. Using 
this factor we can quantify states shown in Fig. 1 , noting 
first a general qualitative agreement with experimental 
findings Q. The case V = 9E r seems almost fully su- 
perfluid (with, however, strongly subpoissonian statistics 
[24J at each site), the case V — 13E r shows first traces 
of insulator phase (integer occupation of sites with van- 
ishing variance, the transition is completed for significant 
fraction of sites &tV — 16E r while for the deepest lattice 
V = 22E r SF fraction is restricted to very narrow regions 
separating different integer occupations. A careful reader 
may notice small indents visible in lines in panels (c) and 
(d) of Fig. [3J They are due to anisotropy induced by a 
cubic lattice in an isotropic harmonic trap. The data, 
presented necessarily as a function of the radial distance 
contain occupations along the axes, the diagonals, and 
other not equivalent directions in the cubic lattice. The 
smallness of the indents indicates that the symmetry of 
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Figure 2: The rescaled mean field momentum distribution (in arbitrary units) in the first Brillouin zone in the (0,0,1) direction 
derived from f° r data presented in Fig.0 The distributions are to be compared with the corresponding "exact" distributions 
in Fig. 2 of [lij]. Observe that mean field result differs from the exact distribution for the case (d) only - see discussion in the 
text. 



the trap is dominant. 



III. TIME-DEPENDENT MEAN FIELD 
DYNAMICS 

The results obtained for mean field ground state in 
realistic situations, shown in the previous Section, seem 
quite encouraging. Yet, in themselves they can say lit- 
tle about the dynamics of the system when the lattice 
depth V is varied. In an attempt to address this impor- 
tant issue we shall use a time-dependent version of the 
mean field approximation. To this end we employ the 
time-dependent variational principle 8] looking for the 
minimum of 

<G(t)\itl^- t -H(t)+nN\G(t)>, (5) 

with H(t) being now the time dependent Hamiltonian. 
The time dependence is implicit via the dependence of 
the BH Hamiltonian H on U, J and Wi, that in turn de- 
pend on V. The chemical potential /x becomes also time 
dependent when system parameters are varied. \G(t) >, 
the variational wavefunction, is assumed in the standard 
Gutzwiller-type form J2J , with /„ (£) now being time- 
dependent. The very same approach has been successfully 
applied recently to the formation of molecules IE Q , the 
treatment of the disordered optical lattices [l(j as well 
as for determining the phase diagram in Bose-Fermi mix- 
tures J13l |. 

The minimization of JSJ yields the set of first order 
differential equations for fn(t) : 



dt In 



u 



J 



n(n- l) + n(Wi - n) 



^yfr+if™ x + $i^t:U , (6) 



(i) 



where $i = Yl<j> < ^(*)l°il^(*) > ( tne sum > as indi- 
cated by subscript in brackets is over the nearest neigh- 
bors only). The nice feature of the evolution resulting 
from equations © is that the average number of parti- 
cles N —< N > is an exact constant of the motion p|. 
Naturally when the parameters of the BH model, e.g. U 
and J, change the chemical potential corresponding to 
the mean field solution with a given number of particles 
N also changes. The dynamics of /i cannot be obtained 
from 10 only. One can find it, however, following the 
evolution of two states \G\ > and |G?2 > with slightly dif- 
ferent average number of particles A^ =< G?2 |7V|G?2 >= 
Ni + SN =< Gi\N\Gi > +6N. The chemical poten- 
tial at given t may be then approximated by fi(t) = (< 
G 2 {t)\H{t)\G 2 {t) > - < GjmH(t)\Gi(t) >)/SN and 
adjusted at each time step 13] . This is the approach 
used in the numerical results presented below. 

Since we want to follow as closely as possible the exper- 
iment [J] let us recall its main features. The experiment 
has three stages after loading the harmonic trap with Rb 
condensate - compare Fig. QJ Firstly, the optical lattice 
depth V(t) was increased in 80 ms (using exponential 
ramp with time constant r = 20 ms) from the initial 
zero value (when the harmonic trap was present only) to 
V m ax = 22E r , where E r is the recoil energy of Rb atoms. 
The sample was then held for 20ms at V max . Finally V(t) 
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Figure 3: Solid lines represent mean field atom density distribution (on-site filling factor - n) as a function of the distance 
from the center of the trap (measured in units of the lattice constant). Dashed lines represent the corresponding variances at 
different sites 2of = 2(< nf > — < n,: > 2 ). The total number of atoms and the SF factor 7sf (see text) for each plot are (a): 

V = 9E r , N = 99771, 7sf = 0.95; (b): V = 13E r , N = 99502, j SF = 0.40; (c): V = lQE r , N = 95408, 7sf = 0.11; (d): 

V — 22E r , N — 94172, 7sf = 0.01. Small indents seen in the lines in panels (c) and (d) are due to anisotropy induced by a 
cubic lattice in an isotropic (harmonic) trap, for further details see text. 



was decreased with the linear ramp to Vf = 9E r with dif- 
ferent speed. At any stage the experiment could be inter- 
rupted by rapidly switching off all laser beams building 
up the lattice as well as the magnetic trap. The freely 
expanding atomic cloud, after some delay, was recorded 
by a destructive absorption imaging, yielding the signal 
which reflects the momentum distribution (141 |lfj . Since 
the absorption images are taken along two orthogonal 
axes the quantity measured is in fact the integrated mo- 
mentum distribution [14j : 



N(k x ,k y ) ex / dk z iik 



(7) 



For clouds released from low optical lattices when tun- 
neling dominates and the superfluid behavior is expected 
the signal reflects Bragg peaks due to interferences of 
atoms coming from different lattice sites. At increased 
lattice depths above 13E r the interference maxima be- 
come immersed in an incoherent background disappear- 
ing practically at 20E r . This behavior was associated 
with the quantum phase transition from SF to MI phase 
U . Most interestingly the coherence of the sample may 



be rapidly recovered when the lattice depth is decreased 
(third stage of the experiment) as measured by the width 
of the central interference peak which decreases almost 
to its original value at V — 9E r in about 4ms. 

In the former applications of the time-dependent mean 
field approach [3, HJ EH 113 the time-dependence of sys- 
tem's parameters was assumed to be sufficiently slow to 
assure adiabaticity. In effect the mean field ground state 
has been followed by applying the time-dependent equa- 
tions for /„ (i)'s ©. Here we have a similar situation 
since it is claimed 4] that the changes in time of V(t) are 
made sufficiently slow to keep the system in the many- 
body ground state. Having the mean field ground states 
for different V values we can (within the mean approxi- 
mation) test this adiabaticity assumption. 

Looking again at the time profile depicted in Fig.0|onc 
may notice that the ramp used in the experiment leads 
indeed to a very slow increase of V(t) initially, however 
changes of V(i) become relatively rapid about and above 
V = 9-EV, i.e. in the region where the transition from 
SF to MI is supposed to take place. Taking as the initial 
state the mean field ground state at V = 9E r , our Simula- 
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Figure 4: The experimental time profile (solid line) of the 
lattice potential depth V measured in the recoil energy E r 
units. The initial exponential increase with a time constant 
t = 20 ms is followed by a flat part and a subsequent decrease 
to Vf = 9£V with a linear ramp with varying slope. Dashed 
line corresponds to the exponential time constant r — 40 ms, 
dash-dotted r = 80 ms. Thin dotted lines indicate particu- 
larly interesting values of time and V - see text for discussion. 



tions show that to assure adiabaticity a small change of V 
on the superfluid side (say from V = 9E r to V = 9.1E r ) 
requires about 20ms (one needs 40 ms for a loop from 
V = 9E r to V = 9.1E r and back to keep the overlap on 
the initial state of the order of 99 percent). That strongly 
indicates that a much longer time is needed to traverse 
adiabatically the whole interesting region from V = 9E r 
to V = 22E r . And that change is realized in about 20ms 
in the experiment. 

To test the adiabatic issue further we shall concentrate 
in the following on the regime above V — 9E r contain- 
ing the quantum phase transition. Starting again from 
the Gutzwiller mean field ground state at V = 9E r we 
simulate the time evolution up to V — 22E r (with ex- 
perimental time profile). We may compare the dynam- 
ically obtained wavefunction plotted in Fig. |SJa) with 
the mean field ground state at V = 22E r (bottom right 
panel (d) in the figure). While the ground state has an 
insulator character almost everywhere in the trap with 
7sf = 0.01, the dynamically evolved wavefunction, by 
comparison, seems to reflect an excited wavepacket and 
it has rather small regions where the occupation of sites 
is close to integer with vanishing number variance. The 
corresponding 75 p = 0.12 confirms the presence of a rel- 
atively large superfluid region. 

To show that the effect is really due to the too fast 
increase of the lattice depth we have modified the ex- 
perimental time profile slightly, by changing the expo- 
nential constant r from 20 ms to 40ms (or 80 ms). That 



makes the initial rise of the laser intensity (and the lattice 
depth) more uniform in time - compare Fig. ^ Observe 
that while the full duration of the first stage remains the 
same, the interval of time spend on the increase of the 
lattice depth from V = 9E r to 22E r increases from below 
20 ms (experimental profile), to about 30 ms (for r = 40 
ms) or about 37 ms (for r = 80 ms) . Starting again from 
the mean field ground state at V = 9E r we obtain the 
atom density profiles shown in Fig.[S{b) and Fig.|3{c), re- 
spectively. Observe that the regions of insulator behavior 
for both cases are much larger than observed previously. 
The corresponding superfluid factors are 7sf = 0.062 for 
t = 40ms and 7sf = 0.048 for t = 80ms. While the 
final distributions still show signs of significant excita- 
tions, the insulator character becomes dominant for (c) 
and also for (b) case. 

Keeping the overall duration of the laser intensity turn- 
on at 80 ms and enlarging the final stage comes at the 
price that the initial rise from V = 0E r to V = 9E r , very 
smooth in the experiment [j] for r = 20 ms, becomes 
sharper for larger r (compare Fig. 0}. Thus larger r 
may lead to some excitation at the initial creation of the 
lattice, not apparent in our simulations since we start 
from the ground state at V = 9E r . Trying to keep the 
total duration of the experiment as short as possible (to 
avoid, for example, the decoherence) one can still imagine 
a slightly more sophisticated pulse rise with say r = 20 
ms initially up to say V = 9E r and further increase with 
a larger time constant, say 40 ms. While the duration of 
the experiment increases by 10% only, the degree of the 
excitation of the final wavepacket becomes much smaller 
and the insulator character much more pronounced. 

In the experiment [J] the atomic density distribution 
is not measured directly. The presence of the Mott insu- 
lator layer has been detected by observing the resonance 
in the excitation spectrum around the interaction energy 
U. Clearly the size of the corresponding peak is related 
to the number of atoms in the insulator layers. Our re- 
sults indicate that an appropriate suggested change in 
the time profile of V(t) should increase the size of insu- 
lator regions and thus enhance the resonant peak in the 
excitation spectrum. 

Let us mention also that qualitatively similar results to 
that depicted in Fig.[5]are obtained starting the evolution 
from different value of V say V = 8E r as long as the ini- 
tial V is sufficiently big but smaller than the beginning of 
the SF-MI transition regime (around V — 13E r ). In the 
experiment the lattice depth is increased from the initial 
zero value but for that case the Bosc-Hubbard Hamilto- 
nian (JIJ is not a good model to describe the cold atoms 
physics. This discrete lattice model is appropriate plUq 
when the atoms have to remain in the lowest vibrational 
state at each site that is for a sufficiently deep lattice. 

It is most interesting to compare the integrated mo- 
mentum distributions J7J) corresponding to the position 
distributions shown in Fig. [S] The cuts along k y = 
are presented in Fig. |SJ In contrast to Fig. [21 we in- 
clude now the Fourier transform of the Wannier function 
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Figure 5: Atom density distribution (on-site filling factor) after a mean field evolution starting from the mean field ground state 
at V = 9E r at the final value of V = 22E r for different V(t) time dependence, (a) corresponds to the exponential time scale of 
20 ms, (b) to 40 ms, (c) to 80 ms. The longer the time scale, the slower the change of V(t) in the considered range as can be 
seen in Fig. [I] Thick dashed lines in (a) to (c) present twice the variance of the on-site occupation showing that the insulator 
regions in (b) and (c) are much larger than in (a). Panel (d) repeats, after Fig. El the mean field ground state at V — 22E r for 
comparison. Observe that indents in the lines for dynamically evolved wavefunctions are much more pronounced than in the 
static (d) case. As before these indents are due to different symmetries of the cubic lattice and sperically harmonic isotropic 
trap. The conflicting symmetries lead to anisotropic excitations in the dynamical situations when the potential parameters are 
varied. Those excitations are responsible for the indents. 



[compare eq. 0J] so the quantity plotted matches the ex- 
periment (for a better quantitative picture, we prefer the 
cut along k y — instead of the three-dimensional color 
plot used oroginally to represent the experimental data) . 
Panel (a) corresponds to the experimental profile. Ob- 
serve that the integrated momentum distribution consists 
of a broad peak in a very nice agreement with the exper- 
iment |4|. Note, however, that for other, "more adia- 
batic" time profiles the narrow central structure emerges 
[compare (b)] becoming quite pronounced both for (c) as 
well as for the mean field ground state momentum dis- 
tribution - depicted in panel (d). The absence of the 
narrow Bragg peak in the experiment seems to be, there- 
fore, not related to the transition to insulator phase as 
suggested in [4j ■ The persistence of the narrow peaks in 
Mott regime has been already noted in [lj| in a model 
static quantum Monte Carlo study on a smaller lattice. 
Result presented in (d) indicates that the conclusion of 
the authors [Ij| concerning the ground state momentum 
distribution extends also to a realistic sample. The dy- 



namical results presented in Fug. |S] suggests that fading 
of the Bragg peaks and the appearance of the broad inte- 
grated momentum distribution can be associated with a 
significant excitation of the sample rather than the tran- 
sition to Mott insulator regime. 

The reliability of the dynamical mean field approach 
may be tested further in an attempt to reproduce 
"restoration of coherence" part of the experiment Q|. 
In the experiment, after reaching V = 22E r the lattice 
depth is kept constant for 20ms and then decreased back 
to V = 9E r with a linear ramp of durationi. It is shown 
that the time needed to restore the narrow interference 
pattern in the integrated momentum distribution is of 
the order of 4ms. This phenomenon is associated with 
restoring the coherence in the sample. 

This interpretation is in contradiction with results al- 
ready presented in Fig. - the existence (or the lack of 
it) of the narrow peaks seems not to be solely related to 
the coherence of the sample but also to its excitation. 
The question, however, remains whether the dynamical 
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Figure 6: Integrated momentum distributions J7J for states shown in Fig. [S] (panels match one another). The cut N(k x ,0) 
along k y — is shown only. Observe that for the time profile of the experiment no narrow peak is visible, in agreement with 
the experiment. For more adiabatic time scale - panel (c), as well as for the ground state - (d) - the narrow peak in the 
momentum distribution is clearly visible. For a further discussion see text. 



mean field calculation can reproduce the results of the 
experiment, i.e. the dependence of the width of the in- 
terference pattern on the time duration of lowering the 
lattice. 

To answer this question we make a simulation (the 
results are shown in Fig. 0), starting from the static 
solution in the SF regime (taken for the convenience 
at V = 9E r again) increasing exponentially the lattice 
height as in the experiment |j], the subsequent delay of 
20 ms at V = 22E r and a linear ramp-down with various 
slopes. Note that the shape of the curve as well as the 
time scale of restoring the narrow interference pattern is 
in quite good agreement with the experiment. While the 
experimental data could be fit with a double exponential 
decay with two time scales, our mean field data are rea- 
sonably reproduced with a single exponential decay with 
time scale r = 1.45 ms. This nicely corresponds with 
the shorter time scale of the experiment (0.94 ms). The 
obtained time scale is also of the order of a typical single 
tunneling time (1/J in appropriate units) to the nearby 
site. Our numerical data fail to reproduce the second ex- 
perimental time scale. This points out to the possibility 
that it can be associated indeed with a long range cor- 
relation between sites. Such a a long range correlation 
should not manifest itself in our mean field simulations 



- the Gutzwiller wavefunction neglects entanglement be- 
tween sites. 

The observed quite good agreement of the obtained 
widths of the momentum distribution with the experi- 
ment J4| seems to be quite a spectacular success of the 
dynamical mean field simulation bearing in mind its sim- 
plicity. The fact that the mean field approach works so 
well may be, in our opinion, attributed to the fact that 
the dynamics takes place in the regime where superfluid 
fraction remains significant. Then the mean atomic field 
$i does not vanish allowing for semiclassical (mean-field) 
description. Our results suggest that the system has 
quite a long memory and remembers that it was originally 
a SF. This fits nicely with the excited wavepacket-like 
character of dynamically obtained wavefunction clearly 
visible in Fig. Ufa). 



IV. CONCLUSIONS 

To summarize, it has been shown that the mean field 
Gutzwiller approximation allows one to simulate a dy- 
namics of inhomogeneous Bose-Hubbard model taking 
into account realistic experimental conditions. The ac- 
curacy of the approximation cannot be controlled which 
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Figure 7: Half-width of the central interference peak for dif- 
ferent ramp down times t r obtained by a lorentzian fit of the 
integrated momentum distribution - compare Q. Filled cir- 
cles are connected by a line to guide the eye. Dashed line is 
a single exponential decay with a time constant r = 1.45 ms. 



is the major drawback of the present approach (a com- 
parison with exact dynamics for small systems will lead 
us nowhere since then the mean field approach is known 
to fail) . On the other hand a comparison with the avail- 



able data seems quite encouraging. Accepting mean field 
predictions we may confirm that indeed the transition 
from superfluid to Mott insulator takes place (albeit in a 
small part of the sample) in the experiment [4J . On the 
other hand the claim that the first stage of the experi- 
ment is performed adiabatically assuring that the system 
remains in its many body ground state (and thus a gen- 
uine textbook quantum phase transition [3j is realized) 
seems questionable. We confirm the suggestion of [lj| 
that fading of narrow peaks in the momentum distribu- 
tion should not be assotiated with the transition to the 
Mott insulator regime. Rather it is a dynamical effect. 

We suggest that optimization of the lattice depth time 
dependence (i.e. laser intensity profile) may help to en- 
large the insulator regions making the transition more 
adiabatic. That may be detected by measuring the size 
of the peak in the excitation spectrum of the system. 

Lastly, let us mention, that a very recent preprint J25J 
reports a study of exact dynamics of the model using the 
method of [20ll21| . However, the results consider at most 
49 atoms in 40 sites of ID lattice. 
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Polish Committee for Scientific Research Grant Quan- 
tum Information and Quantum Engineering, PBZ-MIN- 
008/P03/2003. 
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